Hands-On Exercise 02

Author

Henry Low

Published

September 1, 2024

Modified

September 5, 2024

Data Sources

(saved under ‘data’ folder)
From data.gov.sg:
- Master Plan 2014 Subzone Boundary (Web) (SHP)
- Master Plan 2014 Subzone Boundary (Web) (KML)
- Childcare

Others:
- CostalOutline from SLA

Chapter 4: 1st Order Spatial Point Patterns Analysis Methods

4.1 Setting Up

4.1.1 Loading the R packages

pacman::p_load(sf, raster, spatstat, tmap, tidyverse)

4.1.2 Importing spatial data

# Load childcare dataset
childcare_sf <- st_read("data/child-care-services-geojson.geojson") %>%
  st_transform(crs = 3414) # Project to SVY21
Reading layer `child-care-services-geojson' from data source 
  `C:\Users\Henry\Desktop\SMU Masters\2024-2025 T1\Geospatial Analytics & Applications\Project\GeospatialWebsite\Hands-On_Ex\Hands-On_Ex_02\data\child-care-services-geojson.geojson' 
  using driver `GeoJSON'
Simple feature collection with 1545 features and 2 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 103.6824 ymin: 1.248403 xmax: 103.9897 ymax: 1.462134
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
# Load costaloutline dataset
sg_sf <- st_read(dsn = "data/CostalOutline", layer="CostalOutline") %>%
  st_set_crs(., 3414) # Change CRS attribute to 3414
Reading layer `CostalOutline' from data source 
  `C:\Users\Henry\Desktop\SMU Masters\2024-2025 T1\Geospatial Analytics & Applications\Project\GeospatialWebsite\Hands-On_Ex\Hands-On_Ex_02\data\CostalOutline' 
  using driver `ESRI Shapefile'
Simple feature collection with 60 features and 4 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 2663.926 ymin: 16357.98 xmax: 56047.79 ymax: 50244.03
Projected CRS: SVY21
Warning: st_crs<- : replacing crs does not reproject data; use st_transform for
that
# Load MP14 Subzone dataset
mpsz_sf <- st_read(dsn = "data/MasterPlan2014SubzoneBoundaryWebSHP", layer = "MP14_SUBZONE_WEB_PL") %>%
  st_set_crs(., 3414) # Change CRS attribute to 3414
Reading layer `MP14_SUBZONE_WEB_PL' from data source 
  `C:\Users\Henry\Desktop\SMU Masters\2024-2025 T1\Geospatial Analytics & Applications\Project\GeospatialWebsite\Hands-On_Ex\Hands-On_Ex_02\data\MasterPlan2014SubzoneBoundaryWebSHP' 
  using driver `ESRI Shapefile'
Simple feature collection with 323 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21
Warning: st_crs<- : replacing crs does not reproject data; use st_transform for
that

4.2 Exploratory Plots

# Map MP14 Sf with Childcare Sf 
tm_shape(mpsz_sf) + 
  tm_polygons() + 
tm_shape(childcare_sf) + 
  tm_bubbles(size = 0.01, col = "black")

# Dynamic map with view tmap_mode (requires leaflet package)
tmap_mode('view')
tmap mode set to interactive viewing
tm_shape(childcare_sf)+
  tm_dots()
# Set tmap_mode back to plot
tmap_mode('plot')
tmap mode set to plotting

4.3 Geospatial Data Wrangling

4.3.1 Convert Sf dataframes into Spatial* class

# Convert sf dataframe to spatial class
childcare <- as_Spatial(childcare_sf)
mpsz <- as_Spatial(mpsz_sf)
sg <- as_Spatial(sg_sf)
# Information of childcare spatial class
print(childcare)
class       : SpatialPointsDataFrame 
features    : 1545 
extent      : 11203.01, 45404.24, 25667.6, 49300.88  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
variables   : 2
names       :    Name,                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                           Description 
min values  :   kml_1, <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>ADDRESSBLOCKHOUSENUMBER</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSBUILDINGNAME</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSPOSTALCODE</th> <td>018989</td> </tr><tr bgcolor=""> <th>ADDRESSSTREETNAME</th> <td>1, MARINA BOULEVARD, #B1 - 01, ONE MARINA BOULEVARD, SINGAPORE 018989</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSTYPE</th> <td></td> </tr><tr bgcolor=""> <th>DESCRIPTION</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>HYPERLINK</th> <td></td> </tr><tr bgcolor=""> <th>LANDXADDRESSPOINT</th> <td>0</td> </tr><tr bgcolor="#E3E3F3"> <th>LANDYADDRESSPOINT</th> <td>0</td> </tr><tr bgcolor=""> <th>NAME</th> <td>THE LITTLE SKOOL-HOUSE INTERNATIONAL PTE. LTD.</td> </tr><tr bgcolor="#E3E3F3"> <th>PHOTOURL</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSFLOORNUMBER</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>08F73931F4A691F4</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20200826094036</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSUNITNUMBER</th> <td></td> </tr></table></center> 
max values  : kml_999,                  <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>ADDRESSBLOCKHOUSENUMBER</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSBUILDINGNAME</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSPOSTALCODE</th> <td>829646</td> </tr><tr bgcolor=""> <th>ADDRESSSTREETNAME</th> <td>200, PONGGOL SEVENTEENTH AVENUE, SINGAPORE 829646</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSTYPE</th> <td></td> </tr><tr bgcolor=""> <th>DESCRIPTION</th> <td>Child Care Services</td> </tr><tr bgcolor="#E3E3F3"> <th>HYPERLINK</th> <td></td> </tr><tr bgcolor=""> <th>LANDXADDRESSPOINT</th> <td>0</td> </tr><tr bgcolor="#E3E3F3"> <th>LANDYADDRESSPOINT</th> <td>0</td> </tr><tr bgcolor=""> <th>NAME</th> <td>RAFFLES KIDZ @ PUNGGOL PTE LTD</td> </tr><tr bgcolor="#E3E3F3"> <th>PHOTOURL</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSFLOORNUMBER</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>379D017BF244B0FA</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20200826094036</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSUNITNUMBER</th> <td></td> </tr></table></center> 
# Information of MP Subzone spatial class
print(mpsz)
class       : SpatialPolygonsDataFrame 
features    : 323 
extent      : 2667.538, 56396.44, 15748.72, 50256.33  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
variables   : 15
names       : OBJECTID, SUBZONE_NO, SUBZONE_N, SUBZONE_C, CA_IND, PLN_AREA_N, PLN_AREA_C,       REGION_N, REGION_C,          INC_CRC, FMEL_UPD_D,     X_ADDR,     Y_ADDR,    SHAPE_Leng,    SHAPE_Area 
min values  :        1,          1, ADMIRALTY,    AMSZ01,      N, ANG MO KIO,         AM, CENTRAL REGION,       CR, 00F5E30B5C9B7AD8,      16409,  5092.8949,  19579.069, 871.554887798, 39437.9352703 
max values  :      323,         17,    YUNNAN,    YSSZ09,      Y,     YISHUN,         YS,    WEST REGION,       WR, FFCCF172717C2EAF,      16409, 50424.7923, 49552.7904, 68083.9364708,  69748298.792 
# Information of costaloutline spatial class
print(sg)
class       : SpatialPolygonsDataFrame 
features    : 60 
extent      : 2663.926, 56047.79, 16357.98, 50244.03  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
variables   : 4
names       : GDO_GID, MSLINK, MAPID,              COSTAL_NAM 
min values  :       1,      1,     0,             ISLAND LINK 
max values  :      60,     67,     0, SINGAPORE - MAIN ISLAND 

4.3.2 Convert Spatial* class into generic sp format

# Information of costaloutline spatial class
childcare_sp <- as(childcare, "SpatialPoints")
sg_sp <- as(sg, "SpatialPolygons")
# Information of childcare sp object
print(childcare_sp)
class       : SpatialPoints 
features    : 1545 
extent      : 11203.01, 45404.24, 25667.6, 49300.88  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
# Information of costaloutline sp object
print(sg_sp)
class       : SpatialPolygons 
features    : 60 
extent      : 2663.926, 56047.79, 16357.98, 50244.03  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 

Spatial* classes are specific classes for handling points, lines, polygons, grids, etc., each with dedicated methods and structures.
A generic sp object refers to any object created using the sp package, including all Spatial* classes.

4.3.3 Convert generic sp format into spatstat’s ppp format

# Convert and print information of childcare ppp format
childcare_ppp <- as.ppp(childcare_sf)
Warning in as.ppp.sf(childcare_sf): only first attribute column is used for
marks
childcare_ppp
Marked planar point pattern: 1545 points
marks are of storage type  'character'
window: rectangle = [11203.01, 45404.24] x [25667.6, 49300.88] units
# Plot childcare_ppp
plot(childcare_ppp)
Warning in default.charmap(ntypes, chars): Too many types to display every type
as a different character
Warning: Only 10 out of 1545 symbols are shown in the symbol map

# Check summary statistics of childcare_ppp
summary(childcare_ppp)
Marked planar point pattern:  1545 points
Average intensity 1.91145e-06 points per square unit

Coordinates are given to 11 decimal places

marks are of type 'character'
Summary:
   Length     Class      Mode 
     1545 character character 

Window: rectangle = [11203.01, 45404.24] x [25667.6, 49300.88] units
                    (34200 x 23630 units)
Window area = 808287000 square units

4.3.4 Check for duplicates

# Check if childcare_ppp contains duplicates 
any(duplicated(childcare_ppp))
[1] FALSE
# Check number of duplicated point events
sum(multiplicity(childcare_ppp) > 1)
[1] 0
# Check number of duplicated point events
tmap_mode('view')
tmap mode set to interactive viewing
tm_shape(childcare) +
  tm_dots(alpha=0.4, 
          size=0.05)
# Set tmap_mode back to plot
tmap_mode('plot')
tmap mode set to plotting
# Use jittering to remove duplicates (if any)
childcare_ppp_jit <- rjitter(childcare_ppp, 
                             retry=TRUE, 
                             nsim=1, 
                             drop=TRUE)

# Check if contains any duplicated points
any(duplicated(childcare_ppp_jit))
[1] FALSE

4.3.5 Owin Object

# Create owin object from sf object
sg_owin <- as.owin(sg_sf)
# Plot owin output object
plot(sg_owin)

# Show summary statistics of owin object
summary(sg_owin)
Window: polygonal boundary
50 separate polygons (1 hole)
                 vertices         area relative.area
polygon 1 (hole)       30     -7081.18     -9.76e-06
polygon 2              55     82537.90      1.14e-04
polygon 3              90    415092.00      5.72e-04
polygon 4              49     16698.60      2.30e-05
polygon 5              38     24249.20      3.34e-05
polygon 6             976  23344700.00      3.22e-02
polygon 7             721   1927950.00      2.66e-03
polygon 8            1992   9992170.00      1.38e-02
polygon 9             330   1118960.00      1.54e-03
polygon 10            175    925904.00      1.28e-03
polygon 11            115    928394.00      1.28e-03
polygon 12             24      6352.39      8.76e-06
polygon 13            190    202489.00      2.79e-04
polygon 14             37     10170.50      1.40e-05
polygon 15             25     16622.70      2.29e-05
polygon 16             10      2145.07      2.96e-06
polygon 17             66     16184.10      2.23e-05
polygon 18           5195 636837000.00      8.78e-01
polygon 19             76    312332.00      4.31e-04
polygon 20            627  31891300.00      4.40e-02
polygon 21             20     32842.00      4.53e-05
polygon 22             42     55831.70      7.70e-05
polygon 23             67   1313540.00      1.81e-03
polygon 24            734   4690930.00      6.47e-03
polygon 25             16      3194.60      4.40e-06
polygon 26             15      4872.96      6.72e-06
polygon 27             15      4464.20      6.15e-06
polygon 28             14      5466.74      7.54e-06
polygon 29             37      5261.94      7.25e-06
polygon 30            111    662927.00      9.14e-04
polygon 31             69     56313.40      7.76e-05
polygon 32            143    145139.00      2.00e-04
polygon 33            397   2488210.00      3.43e-03
polygon 34             90    115991.00      1.60e-04
polygon 35             98     62682.90      8.64e-05
polygon 36            165    338736.00      4.67e-04
polygon 37            130     94046.50      1.30e-04
polygon 38             93    430642.00      5.94e-04
polygon 39             16      2010.46      2.77e-06
polygon 40            415   3253840.00      4.49e-03
polygon 41             30     10838.20      1.49e-05
polygon 42             53     34400.30      4.74e-05
polygon 43             26      8347.58      1.15e-05
polygon 44             74     58223.40      8.03e-05
polygon 45            327   2169210.00      2.99e-03
polygon 46            177    467446.00      6.44e-04
polygon 47             46    699702.00      9.65e-04
polygon 48              6     16841.00      2.32e-05
polygon 49             13     70087.30      9.66e-05
polygon 50              4      9459.63      1.30e-05
enclosing rectangle: [2663.93, 56047.79] x [16357.98, 50244.03] units
                     (53380 x 33890 units)
Window area = 725376000 square units
Fraction of frame area: 0.401
# Combine point events object and owin object
childcareSG_ppp = childcare_ppp[sg_owin]
# Show summary statistics of combined object
summary(childcareSG_ppp)
Marked planar point pattern:  1545 points
Average intensity 2.129929e-06 points per square unit

Coordinates are given to 11 decimal places

marks are of type 'character'
Summary:
   Length     Class      Mode 
     1545 character character 

Window: polygonal boundary
50 separate polygons (1 hole)
                 vertices         area relative.area
polygon 1 (hole)       30     -7081.18     -9.76e-06
polygon 2              55     82537.90      1.14e-04
polygon 3              90    415092.00      5.72e-04
polygon 4              49     16698.60      2.30e-05
polygon 5              38     24249.20      3.34e-05
polygon 6             976  23344700.00      3.22e-02
polygon 7             721   1927950.00      2.66e-03
polygon 8            1992   9992170.00      1.38e-02
polygon 9             330   1118960.00      1.54e-03
polygon 10            175    925904.00      1.28e-03
polygon 11            115    928394.00      1.28e-03
polygon 12             24      6352.39      8.76e-06
polygon 13            190    202489.00      2.79e-04
polygon 14             37     10170.50      1.40e-05
polygon 15             25     16622.70      2.29e-05
polygon 16             10      2145.07      2.96e-06
polygon 17             66     16184.10      2.23e-05
polygon 18           5195 636837000.00      8.78e-01
polygon 19             76    312332.00      4.31e-04
polygon 20            627  31891300.00      4.40e-02
polygon 21             20     32842.00      4.53e-05
polygon 22             42     55831.70      7.70e-05
polygon 23             67   1313540.00      1.81e-03
polygon 24            734   4690930.00      6.47e-03
polygon 25             16      3194.60      4.40e-06
polygon 26             15      4872.96      6.72e-06
polygon 27             15      4464.20      6.15e-06
polygon 28             14      5466.74      7.54e-06
polygon 29             37      5261.94      7.25e-06
polygon 30            111    662927.00      9.14e-04
polygon 31             69     56313.40      7.76e-05
polygon 32            143    145139.00      2.00e-04
polygon 33            397   2488210.00      3.43e-03
polygon 34             90    115991.00      1.60e-04
polygon 35             98     62682.90      8.64e-05
polygon 36            165    338736.00      4.67e-04
polygon 37            130     94046.50      1.30e-04
polygon 38             93    430642.00      5.94e-04
polygon 39             16      2010.46      2.77e-06
polygon 40            415   3253840.00      4.49e-03
polygon 41             30     10838.20      1.49e-05
polygon 42             53     34400.30      4.74e-05
polygon 43             26      8347.58      1.15e-05
polygon 44             74     58223.40      8.03e-05
polygon 45            327   2169210.00      2.99e-03
polygon 46            177    467446.00      6.44e-04
polygon 47             46    699702.00      9.65e-04
polygon 48              6     16841.00      2.32e-05
polygon 49             13     70087.30      9.66e-05
polygon 50              4      9459.63      1.30e-05
enclosing rectangle: [2663.93, 56047.79] x [16357.98, 50244.03] units
                     (53380 x 33890 units)
Window area = 725376000 square units
Fraction of frame area: 0.401
# Plot newly derived childcareSG_ppp
plot(childcareSG_ppp)
Warning in default.charmap(ntypes, chars): Too many types to display every type
as a different character
Warning: Only 10 out of 1545 symbols are shown in the symbol map

4.4 First-order Spatial Point Patterns Analysis

# Computing kernel density estimation
kde_childcareSG_bw <- density(childcareSG_ppp,
                              sigma=bw.diggle,
                              edge=TRUE,
                            kernel="gaussian")
# Plot derived kernel density
plot(kde_childcareSG_bw)

# Retrieve bandwidth used to compute kde layer
bw <- bw.diggle(childcareSG_ppp)
bw
   sigma 
298.4095 
# Use rescale.ppp to convert unit of measurement from meter to kilometer
childcareSG_ppp.km <- rescale.ppp(childcareSG_ppp, 1000, "km")
# Re-run density with rescaled dataset
kde_childcareSG.bw <- density(childcareSG_ppp.km, sigma=bw.diggle, edge=TRUE, kernel="gaussian")
# Plot updated density
plot(kde_childcareSG.bw)

# Check bandwidth return by other calculation methods
bw.CvL(childcareSG_ppp.km)
   sigma 
4.543278 
cat("------------------------------------------------------------------\n")
------------------------------------------------------------------
bw.scott(childcareSG_ppp.km)
 sigma.x  sigma.y 
2.224898 1.450966 
cat("------------------------------------------------------------------\n")
------------------------------------------------------------------
bw.ppl(childcareSG_ppp.km)
    sigma 
0.3897114 
cat("------------------------------------------------------------------\n")
------------------------------------------------------------------
bw.diggle(childcareSG_ppp.km)
    sigma 
0.2984095 
# Compare output using bw.diggle and bw.ppl methods
kde_childcareSG.ppl <- density(childcareSG_ppp.km, 
                               sigma=bw.ppl, 
                               edge=TRUE,
                               kernel="gaussian")
par(mfrow=c(1,2), mai=c(0, 0.1, 0.1, 0.3),  oma=c(0, 0, 0, 0), cex.main = 0.8, cex.axis = 0.6)
plot(kde_childcareSG.bw, main = "bw.diggle")
plot(kde_childcareSG.ppl, main = "bw.ppl")

# Compare bw.ppl output for other kernel methods  
par(mfrow=c(2,2), mai=c(0, 0.1, 0.1, 0.3),  oma=c(0, 0, 0, 0), cex.main = 0.8, cex.axis = 0.6)
plot(density(childcareSG_ppp.km, 
             sigma=bw.ppl, 
             edge=TRUE, 
             kernel="gaussian"), 
     main="Gaussian")
plot(density(childcareSG_ppp.km, 
             sigma=bw.ppl, 
             edge=TRUE, 
             kernel="epanechnikov"), 
     main="Epanechnikov")
Warning in density.ppp(childcareSG_ppp.km, sigma = bw.ppl, edge = TRUE, :
Bandwidth selection will be based on Gaussian kernel
plot(density(childcareSG_ppp.km, 
             sigma=bw.ppl, 
             edge=TRUE, 
             kernel="quartic"), 
     main="Quartic")
Warning in density.ppp(childcareSG_ppp.km, sigma = bw.ppl, edge = TRUE, :
Bandwidth selection will be based on Gaussian kernel
plot(density(childcareSG_ppp.km, 
             sigma=bw.ppl, 
             edge=TRUE, 
             kernel="disc"), 
     main="Disc")
Warning in density.ppp(childcareSG_ppp.km, sigma = bw.ppl, edge = TRUE, :
Bandwidth selection will be based on Gaussian kernel

4.4.1 Fixed and Adaptive KDE

# Compute kde layer with a fixed bandwidth
kde_childcareSG_600 <- density(childcareSG_ppp.km, sigma=0.6, edge=TRUE, kernel="gaussian")
plot(kde_childcareSG_600)

# Compute kde layer with an adaptive bandwidth
kde_childcareSG_adaptive <- adaptive.density(childcareSG_ppp.km, method="kernel")
plot(kde_childcareSG_adaptive)

# Compare both fixed and adaptive KDE output
par(mfrow=c(1,2), mai=c(0, 0.1, 0.1, 0.3),  oma=c(0, 0, 0, 0), cex.main = 0.8, cex.axis = 0.6)
plot(kde_childcareSG.bw, main = "Fixed bandwidth")
plot(kde_childcareSG_adaptive, main = "Adaptive bandwidth")

# # Convert KDE output into grid object
gridded_kde_childcareSG_bw <- maptools::as.SpatialGridDataFrame.im(kde_childcareSG.bw) # Use maptools
Please note that 'maptools' will be retired during October 2023,
plan transition at your earliest convenience (see
https://r-spatial.org/r/2023/05/15/evolution4.html and earlier blogs
for guidance);some functionality will be moved to 'sp'.
 Checking rgeos availability: FALSE
spplot(gridded_kde_childcareSG_bw)

# Convert kde_childcareSG.bw to raster
kde_childcareSG_bw_raster <- raster(kde_childcareSG.bw)
kde_childcareSG_bw_raster
class      : RasterLayer 
dimensions : 128, 128, 16384  (nrow, ncol, ncell)
resolution : 0.4170614, 0.2647348  (x, y)
extent     : 2.663926, 56.04779, 16.35798, 50.24403  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : memory
names      : layer 
values     : -8.476185e-15, 28.51831  (min, max)
# Assign CRS information to Rasterlayer
projection(kde_childcareSG_bw_raster) <- CRS("+init=EPSG:3414")
kde_childcareSG_bw_raster
class      : RasterLayer 
dimensions : 128, 128, 16384  (nrow, ncol, ncell)
resolution : 0.4170614, 0.2647348  (x, y)
extent     : 2.663926, 56.04779, 16.35798, 50.24403  (xmin, xmax, ymin, ymax)
crs        : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +units=m +no_defs 
source     : memory
names      : layer 
values     : -8.476185e-15, 28.51831  (min, max)
# Visualize raster with tmap
tm_shape(kde_childcareSG_bw_raster) + 
  tm_raster("layer", palette = "viridis") +
  tm_layout(legend.position = c("right", "bottom"), frame = FALSE)

4.4.1.1 Comparing Spatial Point Patterns with KDE

# Extract study area
pg <- mpsz_sf %>%
  filter(PLN_AREA_N == "PUNGGOL")
tm <- mpsz_sf %>%
  filter(PLN_AREA_N == "TAMPINES")
ck <- mpsz_sf %>%
  filter(PLN_AREA_N == "CHOA CHU KANG")
jw <- mpsz_sf %>%
  filter(PLN_AREA_N == "JURONG WEST")

# Plot Ponggol
par(mfrow=c(2,2), cex.main=0.8)
plot(pg, main = "Ponggol")
Warning: plotting the first 9 out of 15 attributes; use max.plot = 15 to plot
all

# Extract Tampines
par(mfrow=c(2,2), cex.main=0.8)
plot(tm, main = "Tampines")
Warning: plotting the first 9 out of 15 attributes; use max.plot = 15 to plot
all

# Extract Choa Chu Kang
par(mfrow=c(2,2), cex.main=0.8)
plot(ck, main = "Choa Chu Kang")
Warning: plotting the first 10 out of 15 attributes; use max.plot = 15 to plot
all

# Extract Jurong West
par(mfrow=c(2,2), cex.main=0.8)
plot(jw, main = "Jurong West")
Warning: plotting the first 9 out of 15 attributes; use max.plot = 15 to plot
all

## Using 
# Convert sf objects into owin objects
pg_owin = as.owin(pg)
tm_owin = as.owin(tm)
ck_owin = as.owin(ck)
jw_owin = as.owin(jw)

# Combine childcare points and study area
childcare_pg_ppp = childcare_ppp_jit[pg_owin]
childcare_tm_ppp = childcare_ppp_jit[tm_owin]
childcare_ck_ppp = childcare_ppp_jit[ck_owin]
childcare_jw_ppp = childcare_ppp_jit[jw_owin]

# Transform unit of measurement from metre to kilometre
childcare_pg_ppp.km = rescale.ppp(childcare_pg_ppp, 1000, "km")
childcare_tm_ppp.km = rescale.ppp(childcare_tm_ppp, 1000, "km")
childcare_ck_ppp.km = rescale.ppp(childcare_ck_ppp, 1000, "km")
childcare_jw_ppp.km = rescale.ppp(childcare_jw_ppp, 1000, "km")

# Plot the 4 study areas and location of childcare centres
par(mfrow=c(2,2),  mai=c(0.2, 0.2, 0.2, 0.2), cex.main=0.8, cex.axis=0.8)

# par(mfrow=c(1,1))  # Reset to single plot layout
# par(cex.main=1, cex.lab=1, cex.axis=1)  # Reset font sizes to default

plot(childcare_pg_ppp.km, main="Punggol", cex = 0.5)
Warning in default.charmap(ntypes, chars): Too many types to display every type
as a different character
Warning: Only 10 out of 61 symbols are shown in the symbol map
plot(childcare_tm_ppp.km, main="Tampines", cex = 0.5)
Warning in default.charmap(ntypes, chars): Too many types to display every type
as a different character
Warning: Only 10 out of 89 symbols are shown in the symbol map
plot(childcare_ck_ppp.km, main="Choa Chu Kang", cex = 0.5)
Warning in default.charmap(ntypes, chars): Too many types to display every type
as a different character
Warning: Only 10 out of 61 symbols are shown in the symbol map
plot(childcare_jw_ppp.km, main="Jurong West", cex = 0.5)
Warning in default.charmap(ntypes, chars): Too many types to display every type
as a different character
Warning: Only 10 out of 88 symbols are shown in the symbol map

# Plot KDE of planning areas with bw.diggle method
par(mfrow=c(2,2),  mai=c(0.2, 0.2, 0.2, 0.2),  oma=c(0, 0, 0, 0), cex.main = 0.8, cex.axis = 0.7)
plot(density(childcare_pg_ppp.km, 
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Punggol")
plot(density(childcare_tm_ppp.km, 
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Tempines")
plot(density(childcare_ck_ppp.km, 
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Choa Chu Kang")
plot(density(childcare_jw_ppp.km, 
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="gaussian"),
     main="JUrong West")

# Plot KDE of planning areas with fixed bandwidth
par(mfrow=c(2,2), mai=c(0.2, 0.2, 0.2, 0.2),  oma=c(0, 0, 0, 0), cex.main = 0.8, cex.axis = 0.7)
plot(density(childcare_ck_ppp.km, 
             sigma=0.25, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Chou Chu Kang")
plot(density(childcare_jw_ppp.km, 
             sigma=0.25, 
             edge=TRUE, 
             kernel="gaussian"),
     main="JUrong West")
plot(density(childcare_pg_ppp.km, 
             sigma=0.25, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Punggol")
plot(density(childcare_tm_ppp.km, 
             sigma=0.25, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Tampines")

4.4.2 Nearest Neighbour Analysis

# Use Clark and Evans Test for full data
clarkevans.test(childcareSG_ppp,
                correction="none",
                clipregion="sg_owin",
                alternative=c("clustered"),
                nsim=99)

    Clark-Evans test
    No edge correction
    Z-test

data:  childcareSG_ppp
R = 0.55631, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)

Conclusion: Points on childcareSG_ppp are more clustered than randomly distributed (Accept H1) at 95% confidence level (p-value < 0.05)

# Use Clark and Evans Test for Choa Chu Kang
clarkevans.test(childcare_ck_ppp,
                correction="none",
                clipregion=NULL,
                alternative=c("two.sided"),
                nsim=999)

    Clark-Evans test
    No edge correction
    Z-test

data:  childcare_ck_ppp
R = 0.92614, p-value = 0.2698
alternative hypothesis: two-sided
# Use Clark and Evans Test for Tampines
clarkevans.test(childcare_tm_ppp,
                correction="none",
                clipregion=NULL,
                alternative=c("two.sided"),
                nsim=999)

    Clark-Evans test
    No edge correction
    Z-test

data:  childcare_tm_ppp
R = 0.81761, p-value = 0.0009954
alternative hypothesis: two-sided

Chapter 5: 2nd Order Spatial Point Patterns Analysis Methods

5.1 Using G-Function

G function measures the distribution of the distances from an arbitrary event to its nearest event

# Compute G-function for Choa Chu Kang
G_CK = Gest(childcare_ck_ppp, correction = "border")
plot(G_CK, xlim=c(0,500))

# Complete Spatial Randomness Test (Monte Carlo)
G_CK.csr <- envelope(childcare_ck_ppp, Gest, nsim = 999)
Generating 999 simulations of CSR  ...
1, 2, 3, ......10.........20.........30.........40.........50.........60..
.......70.........80.........90.........100.........110.........120.........130
.........140.........150.........160.........170.........180.........190........
.200.........210.........220.........230.........240.........250.........260......
...270.........280.........290.........300.........310.........320.........330....
.....340.........350.........360.........370.........380.........390.........400..
.......410.........420.........430.........440.........450.........460.........470
.........480.........490.........500.........510.........520.........530........
.540.........550.........560.........570.........580.........590.........600......
...610.........620.........630.........640.........650.........660.........670....
.....680.........690.........700.........710.........720.........730.........740..
.......750.........760.........770.........780.........790.........800.........810
.........820.........830.........840.........850.........860.........870........
.880.........890.........900.........910.........920.........930.........940......
...950.........960.........970.........980.........990........
999.

Done.
# Plot results
plot(G_CK.csr)

# Compute G-function for Tampines
G_tm = Gest(childcare_tm_ppp, correction = "best")
plot(G_tm)

# Complete Spatial Randomness Test (Monte Carlo)
G_tm.csr <- envelope(childcare_tm_ppp, Gest, correction = "all", nsim = 999)
Generating 999 simulations of CSR  ...
1, 2, 3, ......10 [1:04 remaining] .........20 [35 sec remaining] ...
......30 [26 sec remaining] .........40 [21 sec remaining] .........50 [19 sec remaining] ..
.......60 [17 sec remaining] .........70 [15 sec remaining] .........80 [14 sec remaining] .
........90 [13 sec remaining] .........100 [13 sec remaining] .........110
 [12 sec remaining] .........120 [12 sec remaining] .........130 [12 sec remaining] .........
140 [11 sec remaining] .........150 [11 sec remaining] .........160 [11 sec remaining] ........
.170 [10 sec remaining] .........180 [10 sec remaining] .........190 [10 sec remaining] .......
..200 [10 sec remaining] .........210 [9 sec remaining] .........220 [9 sec remaining] ......
...230 [9 sec remaining] .........240 [9 sec remaining] .........250 [9 sec remaining] .....
....260 [9 sec remaining] .........270 [8 sec remaining] .........280 [8 sec remaining] ....
.....290 [8 sec remaining] .........300 [8 sec remaining] .........310 [8 sec remaining] ...
......320 [8 sec remaining] .........330 [7 sec remaining] .........340 [7 sec remaining] ..
.......350 [7 sec remaining] .........360 [7 sec remaining] .........370 [7 sec remaining] .
........380 [7 sec remaining] .........390 [7 sec remaining] .........400
 [7 sec remaining] .........410 [6 sec remaining] .........420 [6 sec remaining] .........
430 [6 sec remaining] .........440 [6 sec remaining] .........450 [6 sec remaining] ........
.460 [6 sec remaining] .........470 [6 sec remaining] .........480 [6 sec remaining] .......
..490 [6 sec remaining] .........500 [5 sec remaining] .........510 [5 sec remaining] ......
...520 [5 sec remaining] .........530 [5 sec remaining] .........540 [5 sec remaining] .....
....550 [5 sec remaining] .........560 [5 sec remaining] .........570 [5 sec remaining] ....
.....580 [5 sec remaining] .........590 [5 sec remaining] .........600 [4 sec remaining] ...
......610 [4 sec remaining] .........620 [4 sec remaining] .........630 [4 sec remaining] ..
.......640 [4 sec remaining] .........650 [4 sec remaining] .........660 [4 sec remaining] .
........670 [4 sec remaining] .........680 [3 sec remaining] .........690
 [3 sec remaining] .........700 [3 sec remaining] .........710 [3 sec remaining] .........
720 [3 sec remaining] .........730 [3 sec remaining] .........740 [3 sec remaining] ........
.750 [3 sec remaining] .........760 [3 sec remaining] .........770 [2 sec remaining] .......
..780 [2 sec remaining] .........790 [2 sec remaining] .........800 [2 sec remaining] ......
...810 [2 sec remaining] .........820 [2 sec remaining] .........830 [2 sec remaining] .....
....840 [2 sec remaining] .........850 [2 sec remaining] .........860 [1 sec remaining] ....
.....870 [1 sec remaining] .........880 [1 sec remaining] .........890 [1 sec remaining] ...
......900 [1 sec remaining] .........910 [1 sec remaining] .........920 [1 sec remaining] ..
.......930 [1 sec remaining] .........940 [1 sec remaining] .........950 [1 sec remaining] .
........960 [0 sec remaining] .........970 [0 sec remaining] .........980
 [0 sec remaining] .........990 [0 sec remaining] ........
999.

Done.
# Plot results
plot(G_tm.csr)

5.2 Using F-Function

F function estimates the empty space function F(r) or its hazard rate h(r) from a point pattern in a window of arbitrary

# Compute F-function for Choa Chu Kang
F_CK = Fest(childcare_ck_ppp)
plot(F_CK)

# Complete Spatial Randomness Test (Monte Carlo)
F_CK.csr <- envelope(childcare_ck_ppp, Fest, nsim = 999)
Generating 999 simulations of CSR  ...
1, 2, 3, ......10.........20.........30.........40.........50.........60..
.......70.........80.........90.........100.........110.........120.........130
.........140.........150.........160.........170.........180.........190........
.200.........210.........220.........230.........240.........250.........260......
...270.........280.........290.........300.........310.........320.........330....
.....340.........350.........360.........370.........380.........390.........400..
.......410.........420.........430.........440.........450.........460.........470
.........480.........490.........500.........510.........520.........530........
.540.........550.........560.........570.........580.........590.........600......
...610.........620.........630.........640.........650.........660.........670....
.....680.........690.........700.........710.........720.........730.........740..
.......750.........760.........770.........780.........790.........800.........810
.........820.........830.........840.........850.........860.........870........
.880.........890.........900.........910.........920.........930.........940......
...950.........960.........970.........980.........990........
999.

Done.
# Plot results
plot(F_CK.csr)

# Compute F-function for Tampines
F_tm = Fest(childcare_tm_ppp, correction = "best")
plot(F_tm)

# Complete Spatial Randomness Test (Monte Carlo)
F_tm.csr <- envelope(childcare_tm_ppp, Fest, correction = "all", nsim = 999)
Generating 999 simulations of CSR  ...
1, 2, 3, ......10.........20.........30.........40.........50.........60..
.......70.........80.........90.........100.........110.........120.........130
.........140.........150.........160.........170.........180.........190........
.200.........210.........220.........230.........240.........250.........260......
...270.........280.........290.........300.........310.........320.........330....
.....340.........350.........360.........370.........380.........390.........400..
.......410.........420.........430.........440.........450.........460.........470
.........480.........490.........500.........510.........520.........530........
.540.........550.........560.........570.........580.........590.........600......
...610.........620.........630.........640.........650.........660.........670....
.....680.........690.........700.........710.........720.........730.........740..
.......750.........760.........770.........780.........790.........800.........810
.........820.........830.........840.........850.........860.........870........
.880.........890.........900.........910.........920.........930.........940......
...950.........960.........970.........980.........990........
999.

Done.
# Plot results
plot(F_tm.csr)

5.3 Using K-Function

K-function measures the number of events found up to a given distance of any particular event.

# Compute K-function for Choa Chu Kang
K_ck = Kest(childcare_ck_ppp, correction = "Ripley")
plot(K_ck, . -r ~ r, ylab= "K(d)-r", xlab = "d(m)")

# Complete Spatial Randomness Test (Monte Carlo)
K_ck.csr <- envelope(childcare_ck_ppp, Kest, nsim = 99, rank = 1, glocal=TRUE)
Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.
# Plot results
plot(K_ck.csr, . - r ~ r, xlab="d", ylab="K(d)-r")

# Compute K-function for Tampines
K_tm = Kest(childcare_tm_ppp, correction = "Ripley")
plot(K_tm, . -r ~ r, 
     ylab= "K(d)-r", xlab = "d(m)", 
     xlim=c(0,1000))

# Complete Spatial Randomness Test (Monte Carlo)
K_tm.csr <- envelope(childcare_tm_ppp, Kest, nsim = 99, rank = 1, glocal=TRUE)
Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.
# Plot results
plot(K_tm.csr, . - r ~ r, 
     xlab="d", ylab="K(d)-r", xlim=c(0,500))

5.4 Using L-Function

# Compute L-function for Choa Chu Kang
L_ck = Lest(childcare_ck_ppp, correction = "Ripley")
plot(L_ck, . -r ~ r, 
     ylab= "L(d)-r", xlab = "d(m)")

# Complete Spatial Randomness Test (Monte Carlo)
L_ck.csr <- envelope(childcare_ck_ppp, Lest, nsim = 99, rank = 1, glocal=TRUE)
Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.
# Plot results
plot(L_ck.csr, . - r ~ r, xlab="d", ylab="L(d)-r")

# Compute L-function for Tampines
L_tm = Lest(childcare_tm_ppp, correction = "Ripley")
plot(L_tm, . -r ~ r, 
     ylab= "L(d)-r", xlab = "d(m)", 
     xlim=c(0,1000))

# Complete Spatial Randomness Test (Monte Carlo)
L_tm.csr <- envelope(childcare_tm_ppp, Lest, nsim = 99, rank = 1, glocal=TRUE)
Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.
# Plot results
plot(L_tm.csr, . - r ~ r, 
     xlab="d", ylab="L(d)-r", xlim=c(0,500))